knitr::opts_chunk$set(echo = TRUE)
cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
report_function_hash <- digest::digest(summarize_brms)## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE
options(
dplyr.print_max = 100,
brms.backend = 'cmdstan',
brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
brms.file_refit = 'on_change',
#brms.file_refit = 'always',
error = function() {
beepr::beep(sound = 5)
if (shutdown) {
system("shutdown /s /t 180")
quit(save = "no", status = 1)
}
}
, es.use_symbols = TRUE
)
####################### Model parameters #######################
iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000
# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'
################################################################
suffix = as.character(iterations)df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df
df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]Constructing scales reshaping data (4field) centering data within and between
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323
# For indistinguishable Dyads
model_rows_fixed <- c(
'Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw',
'pushing_self_cw',
'pushing_partner_cw',
'day',
'weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'persuasion_self_cb',
'persuasion_partner_cb',
'pressure_self_cb',
'pressure_partner_cb',
'pushing_self_cb',
'pushing_partner_cb',
'weartime_self_cb'
)
model_rows_fixed_ordinal <- c(
model_rows_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rows_fixed[2:length(model_rows_fixed)]
)
model_rows_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(persuasion_self_cw)',
'sd(persuasion_partner_cw)',
'sd(pressure_self_cw)',
'sd(pressure_partner_cw)',
'sd(pushing_self_cw)',
'sd(pushing_partner_cw)',
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rows_random_ordinal <- c(model_rows_random,'disc')# For indistinguishable Dyads
model_rownames_fixed <- c(
"Intercept",
# "-- WITHIN PERSON MAIN EFFECTS --",
"Daily persuasion experienced",
"Daily persuasion utilized (partner's view)", # OR partner received
"Daily pressure experienced",
"Daily pressure utilized (partner's view)",
"Daily pushing experienced",
"Daily pushing utilized (partner's view)",
"Day",
"Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Mean persuasion experienced",
"Mean persuasion utilized (partner's view)",
"Mean pressure experienced",
"Mean pressure utilized (partner's view)",
"Mean pushing experienced",
"Mean pushing utilized (partner's view)",
"Mean weartime"
)
model_rownames_fixed_ordinal <- c(
model_rownames_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed[2:length(model_rownames_fixed)]
)
model_rownames_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
"sd(Daily persuasion experienced)",
"sd(Daily persuasion utilized (partner's view))", # OR partner received
"sd(Daily pressure experienced)",
"sd(Daily pressure utilized (partner's view))",
"sd(Daily pushing experienced)",
"sd(Daily pushing utilized (partner's view))",
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rownames_random_ordinal <- c(model_rownames_random,'disc')rows_to_pack <- list(
"Within-Person Effects" = c(2,9),
"Between-Person Effects" = c(10,16),
"Random Effects" = c(17, 23),
"Additional Parameters" = c(24,24)
)
rows_to_pack_ordinal <- list(
"Intercepts" = c(1,6),
"Within-Person Effects" = c(2+5,9+5),
"Between-Person Effects" = c(10+5,16+5),
"Random Effects" = c(17+5, 23+5),
"Additional Parameters" = c(24+5,24+6)
)HURDLE MODELS
# For indistinguishable Dyads
model_rows_fixed_hu <- c(
'Intercept',
'hu_Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw',
'pushing_self_cw',
'pushing_partner_cw',
'day',
'weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'persuasion_self_cb',
'persuasion_partner_cb',
'pressure_self_cb',
'pressure_partner_cb',
'pushing_self_cb',
'pushing_partner_cb',
'weartime_self_cb',
# HURDLE MODEL
# '-- WITHIN PERSON MAIN EFFECTS --',
'hu_persuasion_self_cw',
'hu_persuasion_partner_cw',
'hu_pressure_self_cw',
'hu_pressure_partner_cw',
'hu_pushing_self_cw',
'hu_pushing_partner_cw',
'hu_day',
'hu_weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'hu_persuasion_self_cb',
'hu_persuasion_partner_cb',
'hu_pressure_self_cb',
'hu_pressure_partner_cb',
'hu_pushing_self_cb',
'hu_pushing_partner_cb',
'hu_weartime_self_cb'
)
model_rows_fixed_hu_ordinal <- c(
model_rows_fixed_hu[1:2],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)
model_rows_random_hu <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(hu_Intercept)',
'sd(persuasion_self_cw)',
'sd(persuasion_partner_cw)',
'sd(pressure_self_cw)',
'sd(pressure_partner_cw)',
'sd(pushing_self_cw)',
'sd(pushing_partner_cw)',
# HURDLE
'sd(hu_persuasion_self_cw)',
'sd(hu_persuasion_partner_cw)',
'sd(hu_pressure_self_cw)',
'sd(hu_pressure_partner_cw)',
'sd(hu_pushing_self_cw)',
'sd(hu_pushing_partner_cw)',
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
"Intercept",
"Hurdle Intercept",
# "-- WITHIN PERSON MAIN EFFECTS --",
"Daily persuasion experienced",
"Daily persuasion utilized (partner's view)", # OR partner received
"Daily pressure experienced",
"Daily pressure utilized (partner's view)",
"Daily pushing experienced",
"Daily pushing utilized (partner's view)",
"Day",
"Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Mean persuasion experienced",
"Mean persuasion utilized (partner's view)",
"Mean pressure experienced",
"Mean pressure utilized (partner's view)",
"Mean pushing experienced",
"Mean pushing utilized (partner's view)",
"Mean weartime",
# HURDLE
# "-- WITHIN PERSON MAIN EFFECTS --",
"Hu Daily persuasion experienced",
"Hu Daily persuasion utilized (partner's view)", # OR partner received
"Hu Daily pressure experienced",
"Hu Daily pressure utilized (partner's view)",
"Hu Daily pushing experienced",
"Hu Daily pushing utilized (partner's view)",
"Hu Day",
"Hu Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Hu Mean persuasion experienced",
"Hu Mean persuasion utilized (partner's view)",
"Hu Mean pressure experienced",
"Hu Mean pressure utilized (partner's view)",
"Hu Mean pushing experienced",
"Hu Mean pushing utilized (partner's view)",
"Hu Mean weartime"
)
model_rownames_fixed_hu_ordinal <- c(
model_rownames_fixed_hu[1:2],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)
model_rownames_random_hu <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(Hurdle Intercept)',
"sd(Daily persuasion experienced)",
"sd(Daily persuasion utilized (partner's view))", # OR partner received
"sd(Daily pressure experienced)",
"sd(Daily pressure utilized (partner's view))",
"sd(Daily pushing experienced)",
"sd(Daily pushing utilized (partner's view))",
# Hurdle
"sd(Hu Daily persuasion experienced)",
"sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
"sd(Hu Daily pressure experienced)",
"sd(Hu Daily pressure utilized (partner's view))",
"sd(Hu Daily pushing experienced)",
"sd(Hu Daily pushing utilized (partner's view))",
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')rows_to_pack_hu <- list(
"Conditional Within-Person Effects" = c(3,10),
"Conditional Between-Person Effects" = c(11,17),
"Hurdle Within-Person Effects" = c(18,25),
"Hurdle Between-Person Effects" = c(26,32),
"Random Effects" = c(33, 46),
"Additional Parameters" = c(47,47)
)
rows_to_pack_hu_ordinal <- list(
"Intercepts" = c(1,7),
"Conditional Within-Person Effects" = c(3+5,10+5),
"Conditional Between-Person Effects" = c(11+5,17+5),
"Hurdle Within-Person Effects" = c(18+5,25+5),
"Hurdle Between-Person Effects" = c(26+5,32+5),
"Random Effects" = c(33+5, 46+5),
"Additional Parameters" = c(47+5,47+6)
)formula <- bf(
pa_sub ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day +
# Random effects
(1 + persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID),
hu = ~ persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day +
# Random effects
(1 + persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
, brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_Nopushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 3736 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -10536.2 166.2
## p_loo 144.3 4.8
## looic 21072.3 332.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.5]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 3735 100.0% 512
## (0.7, 1] (bad) 1 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 8, observations = 3736, p-value = 0.46
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000267666 0.003211991
## sample estimates:
## outlier frequency (expected: 0.00152837259100642 )
## 0.002141328
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
pa_sub,
stats_to_report = c('pd','ROPE'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub_continuous <- summarize_brms(
pa_sub,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_continuous,
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous
columns_to_replace <- c("ROPE", "inside ROPE")
summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <-
summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]
# Print the updated dataframe
summary_pa_sub %>%
print_df(rows_to_pack = rows_to_pack_hu)| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 48.72*** | 3.08 | [42.92, 55.10] | 1.000 | [0.92, 1.08] | 0.000 | 1.002 | 964 | 1749 |
| Hurdle Intercept | 0.82 | 0.14 | [ 0.59, 1.13] | 0.882 | [0.84, 1.20] | 0.433 | 1.005 | 552 | 1254 |
| Conditional Within-Person Effects | |||||||||
| Daily persuasion experienced | 1.03 | 0.03 | [ 0.97, 1.09] | 0.853 | [0.92, 1.08] | 0.965 | 1.001 | 1471 | 2332 |
| Daily persuasion utilized (partner’s view) | 1.03 | 0.02 | [ 0.98, 1.08] | 0.877 | [0.92, 1.08] | 0.987 | 1.002 | 2066 | 2263 |
| Daily pressure experienced | 0.90* | 0.04 | [ 0.81, 0.99] | 0.983 | [0.92, 1.08] | 0.260 | 1.000 | 3707 | 2336 |
| Daily pressure utilized (partner’s view) | 0.94 | 0.04 | [ 0.86, 1.02] | 0.928 | [0.92, 1.08] | 0.634 | 1.002 | 4371 | 3030 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.98 | 0.06 | [ 0.87, 1.11] | 0.623 | [0.92, 1.08] | 0.766 | 1.000 | 5745 | 2727 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||
| Mean persuasion experienced | 1.11 | 0.15 | [ 0.85, 1.47] | 0.775 | [0.92, 1.08] | 0.336 | 1.004 | 817 | 1463 |
| Mean persuasion utilized (partner’s view) | 1.09 | 0.15 | [ 0.83, 1.45] | 0.737 | [0.92, 1.08] | 0.354 | 1.007 | 812 | 1214 |
| Mean pressure experienced | 1.12 | 0.18 | [ 0.81, 1.55] | 0.752 | [0.92, 1.08] | 0.287 | 1.005 | 1015 | 1888 |
| Mean pressure utilized (partner’s view) | 0.91 | 0.16 | [ 0.64, 1.28] | 0.708 | [0.92, 1.08] | 0.305 | 1.003 | 1145 | 1670 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||
| Hu Daily persuasion experienced | 1.61*** | 0.10 | [ 1.43, 1.84] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 1969 | 2771 |
| Hu Daily persuasion utilized (partner’s view) | 1.43*** | 0.08 | [ 1.28, 1.61] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3179 | 3039 |
| Hu Daily pressure experienced | 1.01 | 0.15 | [ 0.76, 1.39] | 0.525 | [0.84, 1.20] | 0.779 | 1.000 | 3515 | 2678 |
| Hu Daily pressure utilized (partner’s view) | 1.74*** | 0.32 | [ 1.26, 2.83] | 1.000 | [0.84, 1.20] | 0.009 | 1.001 | 2286 | 1405 |
| Hu Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.92 | 0.12 | [ 0.71, 1.20] | 0.730 | [0.84, 1.20] | 0.759 | 1.001 | 7057 | 2959 |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||
| Hu Mean persuasion experienced | 1.56 | 0.55 | [ 0.76, 3.18] | 0.889 | [0.84, 1.20] | 0.183 | 1.011 | 579 | 1175 |
| Hu Mean persuasion utilized (partner’s view) | 1.55 | 0.55 | [ 0.76, 3.19] | 0.881 | [0.84, 1.20] | 0.188 | 1.009 | 606 | 1142 |
| Hu Mean pressure experienced | 0.33** | 0.14 | [ 0.15, 0.77] | 0.996 | [0.84, 1.20] | 0.013 | 1.007 | 782 | 1578 |
| Hu Mean pressure utilized (partner’s view) | 0.58 | 0.24 | [ 0.25, 1.33] | 0.911 | [0.84, 1.20] | 0.139 | 1.007 | 803 | 1296 |
| Hu Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.32 | 0.04 | [0.24, 0.42] | NA | NA | NA | 1.008 | 985 | 1765 |
| sd(Hurdle Intercept) | 0.89 | 0.12 | [0.70, 1.18] | NA | NA | NA | 1.001 | 902 | 1729 |
| sd(Daily persuasion experienced) | 0.12 | 0.02 | [0.08, 0.17] | NA | NA | NA | 1.003 | 1813 | 2949 |
| sd(Daily persuasion utilized (partner’s view)) | 0.09 | 0.02 | [0.05, 0.14] | NA | NA | NA | 1.006 | 2186 | 2589 |
| sd(Daily pressure experienced) | 0.08 | 0.07 | [0.00, 0.25] | NA | NA | NA | 1.001 | 1377 | 1876 |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | 0.05 | [0.00, 0.19] | NA | NA | NA | 1.003 | 1648 | 1853 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion experienced) | 0.21 | 0.07 | [0.08, 0.37] | NA | NA | NA | 1.002 | 1341 | 977 |
| sd(Hu Daily persuasion utilized (partner’s view)) | 0.19 | 0.07 | [0.05, 0.34] | NA | NA | NA | 1.001 | 1429 | 1309 |
| sd(Hu Daily pressure experienced) | 0.23 | 0.20 | [0.01, 0.80] | NA | NA | NA | 1.003 | 1209 | 1877 |
| sd(Hu Daily pressure utilized (partner’s view)) | 0.31 | 0.25 | [0.02, 1.01] | NA | NA | NA | 1.004 | 1222 | 1709 |
| sd(Hu Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | 0.69 | 0.01 | [0.67, 0.71] | NA | NA | NA | 1.001 | 5060 | 3213 |
# Plot continuous part of model
variable <- c(
'(Intercept)',
'b_persuasion_self_cw',
'b_persuasion_partner_cw',
'b_pressure_self_cw',
'b_pressure_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()plot(
bayestestR::rope(
pa_sub,
parameter = variable,
range = rope_range_continuous,
verbose = F,
ci = 1
)
) + theme_bw()# Hurdle part of the model
variable <- c(
'b_hu_persuasion_self_cw',
'b_hu_persuasion_partner_cw',
'b_hu_pressure_self_cw',
'b_hu_pressure_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
verbose = FALSE
) + theme_bw()## Possible multicollinearity between b_persuasion_partner_cb and
## b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
## b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
## b_persuasion_partner_cb (r = 0.77), b_hu_persuasion_partner_cb and
## b_hu_persuasion_self_cb (r = 0.82), b_hu_pressure_self_cb and
## b_hu_persuasion_self_cb (r = 0.81), b_hu_pressure_partner_cb and
## b_hu_persuasion_partner_cb (r = 0.81). This might lead to inappropriate
## results. See 'Details' in '?rope'.
conds_eff <- conditional_spaghetti(
pa_sub,
effects = c(
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw'
),
x_label = c(
'Received Persuasion',
'Exerted Persuasion',
'Received Pressure',
'Exerted Pressure'
),
group_var = 'coupleID',
plot_full_range = TRUE,
y_limits = c(0, 100),
y_label = "Same-Day MVPA",
y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
, filter_quantiles = .9995
, font_family = 'Candara'
)## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
$persuasion_self_cw
## Warning: Removed 157 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 173 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.011
$persuasion_partner_cw
## Warning: Removed 143 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 145 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00867
$pressure_self_cw
## Warning: Removed 46 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 55 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.016
$pressure_partner_cw
## Warning: Removed 49 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 85 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0275
Note. This graphic illustrates the relationship between
social control and moderate to vigorous physical activity (MVPA) using a
Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered
within individuals to examine how deviations from their average social
control relate to same-day MVPA. Shaded areas indicate credible
intervals, thick lines show fixed effects, and thin lines represent
random effects, highlighting variability across couples. The plots
display the probability of being active, expected minutes of MVPA when
active, and combined predicted MVPA. The bottom density plot visualizes
the posterior distributions of slope estimates, transformed to represent
multiplicative changes in odds ratios (hurdle component) or expected
values. Medians and 95% credible intervals (2.5th and 97.5th
percentiles) are shown. Effects are significant, when the 95% credible
interval does not overlap 1.
## [1] 5.75 971.25
formula <- bf(
pa_obj ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day + weartime_self_cw + weartime_self_cb +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "Intercept")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_nopushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## Computed from 4000 by 3337 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -18783.3 68.4
## p_loo 74.9 3.4
## looic 37566.6 136.7
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.7]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 27, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.001498352 0.005394067
## sample estimates:
## outlier frequency (expected: 0.00329937069223854 )
## 0.0080911
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summarize_brms(
pa_obj_log,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_log,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 116.98*** | 6.11 | [105.53, 130.15] | 1.000 | [0.94, 1.07] | 0.000 | 1.003 | 583 | 1345 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 1.04* | 0.01 | [ 1.01, 1.07] | 0.989 | [0.94, 1.07] | 0.975 | 1.001 | 2606 | 2982 |
| Daily persuasion utilized (partner’s view) | 1.02 | 0.02 | [ 0.99, 1.06] | 0.936 | [0.94, 1.07] | 0.995 | 1.000 | 3058 | 3179 |
| Daily pressure experienced | 0.95 | 0.03 | [ 0.89, 1.02] | 0.935 | [0.94, 1.07] | 0.678 | 1.000 | 4455 | 3085 |
| Daily pressure utilized (partner’s view) | 0.99 | 0.03 | [ 0.93, 1.06] | 0.627 | [0.94, 1.07] | 0.936 | 1.001 | 4236 | 2961 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.97 | 0.03 | [ 0.91, 1.04] | 0.774 | [0.94, 1.07] | 0.856 | 1.001 | 7048 | 2972 |
| Daily weartime | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 1.001 | 3742 | 2490 |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | 1.11 | 0.14 | [ 0.85, 1.42] | 0.788 | [0.94, 1.07] | 0.279 | 1.007 | 529 | 1039 |
| Mean persuasion utilized (partner’s view) | 1.04 | 0.13 | [ 0.79, 1.34] | 0.611 | [0.94, 1.07] | 0.374 | 1.008 | 538 | 1075 |
| Mean pressure experienced | 0.89 | 0.13 | [ 0.68, 1.20] | 0.782 | [0.94, 1.07] | 0.244 | 1.004 | 669 | 1436 |
| Mean pressure utilized (partner’s view) | 1.02 | 0.14 | [ 0.78, 1.37] | 0.578 | [0.94, 1.07] | 0.362 | 1.005 | 643 | 1395 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | 1.00 | 0.00 | [ 1.00, 1.00] | 0.912 | [0.94, 1.07] | 1.000 | 1.000 | 4781 | 3605 |
| Random Effects | |||||||||
| sd(Intercept) | 0.30 | 0.04 | [0.23, 0.39] | NA | NA | NA | 1.005 | 957 | 1591 |
| sd(Daily persuasion experienced) | 0.05 | 0.01 | [0.03, 0.08] | NA | NA | NA | 1.002 | 2335 | 2032 |
| sd(Daily persuasion utilized (partner’s view)) | 0.06 | 0.02 | [0.03, 0.09] | NA | NA | NA | 1.001 | 1909 | 2207 |
| sd(Daily pressure experienced) | 0.05 | 0.04 | [0.00, 0.14] | NA | NA | NA | 1.001 | 1627 | 2057 |
| sd(Daily pressure utilized (partner’s view)) | 0.04 | 0.03 | [0.00, 0.12] | NA | NA | NA | 1.000 | 2070 | 1630 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | 0.58 | 0.01 | [0.56, 0.59] | NA | NA | NA | 1.004 | 7025 | 2589 |
plot(
bayestestR::p_direction(pa_obj_log),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_persuasion_partner_cb and
## b_persuasion_self_cb (r = 0.88), b_pressure_self_cb and
## b_persuasion_self_cb (r = 0.85), b_pressure_partner_cb and
## b_persuasion_self_cb (r = 0.79), b_pressure_self_cb and
## b_persuasion_partner_cb (r = 0.77), b_pressure_partner_cb and
## b_persuasion_partner_cb (r = 0.87). This might lead to inappropriate
## results. See 'Details' in '?rope'.
## [1] 0 5
formula <- bf(
aff ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_nopushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 3736 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -5189.1 59.2
## p_loo 68.2 3.0
## looic 10378.3 118.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 3735 100.0% 750
## (0.7, 1] (bad) 1 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 32, observations = 3736, p-value =
## 2.24e-11
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.005865832 0.012070325
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.00856531
summarize_brms(
mood_gauss,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)| Est. | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 3.70*** | 0.11 | [ 3.49, 3.90] | 1.000 | [-0.11, 0.11] | 0.000 | 1.006 | 417 | 729 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.01 | 0.02 | [-0.03, 0.05] | 0.663 | [-0.11, 0.11] | 1.000 | 1.001 | 3672 | 2684 |
| Daily persuasion utilized (partner’s view) | 0.04 | 0.02 | [-0.01, 0.08] | 0.938 | [-0.11, 0.11] | 0.999 | 1.000 | 3341 | 3254 |
| Daily pressure experienced | -0.03 | 0.05 | [-0.14, 0.07] | 0.743 | [-0.11, 0.11] | 0.940 | 1.000 | 4889 | 2934 |
| Daily pressure utilized (partner’s view) | 0.01 | 0.05 | [-0.10, 0.11] | 0.559 | [-0.11, 0.11] | 0.965 | 1.001 | 3761 | 2704 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.25*** | 0.06 | [ 0.15, 0.36] | 1.000 | [-0.11, 0.11] | 0.006 | 1.001 | 7106 | 2660 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | 0.42 | 0.25 | [-0.10, 0.91] | 0.942 | [-0.11, 0.11] | 0.099 | 1.010 | 499 | 775 |
| Mean persuasion utilized (partner’s view) | 0.33 | 0.25 | [-0.18, 0.82] | 0.902 | [-0.11, 0.11] | 0.152 | 1.009 | 487 | 558 |
| Mean pressure experienced | -0.36 | 0.26 | [-0.91, 0.19] | 0.903 | [-0.11, 0.11] | 0.139 | 1.007 | 571 | 954 |
| Mean pressure utilized (partner’s view) | -0.29 | 0.27 | [-0.82, 0.27] | 0.851 | [-0.11, 0.11] | 0.180 | 1.007 | 549 | 835 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.61 | 0.07 | [0.48, 0.79] | NA | NA | NA | 1.009 | 812 | 1382 |
| sd(Daily persuasion experienced) | 0.05 | 0.03 | [0.00, 0.11] | NA | NA | NA | 1.005 | 1000 | 1287 |
| sd(Daily persuasion utilized (partner’s view)) | 0.08 | 0.03 | [0.03, 0.14] | NA | NA | NA | 1.003 | 1363 | 1122 |
| sd(Daily pressure experienced) | 0.06 | 0.06 | [0.00, 0.23] | NA | NA | NA | 1.001 | 1944 | 2092 |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | 0.06 | [0.00, 0.25] | NA | NA | NA | 1.002 | 2091 | 2255 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | 0.96 | 0.01 | [0.94, 0.98] | NA | NA | NA | 1.000 | 8235 | 2720 |
plot(
bayestestR::p_direction(mood_gauss),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_pressure_partner_cb and
## b_persuasion_self_cb (r = 0.86), b_pressure_self_cb and
## b_persuasion_partner_cb (r = 0.86), b_pressure_partner_cb and
## b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
## See 'Details' in '?rope'.
## [1] 0 5
df_double$reactance_ordinal <- factor(df_double$reactance,
levels = 0:5,
ordered = TRUE)
formula <- bf(
reactance_ordinal ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = cumulative() # HURDLE_CUMULATIVE
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance_ordinal <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::cumulative(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777
, file = file.path("models_cache_brms", paste0("reactance_ordinal_nopushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 756 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -683.2 31.6
## p_loo 59.6 4.5
## looic 1366.3 63.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 751 99.3% 224
## (0.7, 1] (bad) 5 0.7% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 0.36
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.002645503
## sample estimates:
## outlier frequency (expected: 0.000304232804232804 )
## 0.001322751
summarize_brms(
reactance_ordinal,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_ordinal,
model_rows_random = model_rows_random_ordinal,
model_rownames_fixed = model_rownames_fixed_ordinal,
model_rownames_random = model_rownames_random_ordinal,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack_ordinal)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercepts | |||||||||
| Intercept | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Intercept[1] | 3.48*** | 0.86 | [ 2.16, 5.75] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 3268 | 2969 |
| Intercept[2] | 7.32*** | 1.90 | [ 4.48, 12.47] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3373 | 3061 |
| Intercept[3] | 19.78*** | 5.57 | [ 11.49, 35.14] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3544 | 3388 |
| Intercept[4] | 84.24*** | 29.47 | [ 45.16, 173.12] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 4264 | 3330 |
| Intercept[5] | 2760.30*** | 1715.32 | [929.55, 10049.83] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 5669 | 2968 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.85* | 0.07 | [ 0.72, 0.99] | 0.979 | [0.84, 1.20] | 0.596 | 1.000 | 4158 | 3111 |
| Daily persuasion utilized (partner’s view) | 1.01 | 0.09 | [ 0.84, 1.20] | 0.548 | [0.84, 1.20] | 0.950 | 1.001 | 4250 | 2828 |
| Daily pressure experienced | 1.98** | 0.34 | [ 1.31, 2.76] | 0.999 | [0.84, 1.20] | 0.011 | 1.001 | 2901 | 3115 |
| Daily pressure utilized (partner’s view) | 1.15 | 0.23 | [ 0.72, 1.85] | 0.749 | [0.84, 1.20] | 0.509 | 1.001 | 3088 | 2320 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 1.63 | 0.55 | [ 0.87, 3.10] | 0.929 | [0.84, 1.20] | 0.161 | 1.000 | 6064 | 2790 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | 1.08 | 0.50 | [ 0.42, 2.76] | 0.567 | [0.84, 1.20] | 0.309 | 1.005 | 1276 | 2066 |
| Mean persuasion utilized (partner’s view) | 0.79 | 0.42 | [ 0.29, 2.26] | 0.665 | [0.84, 1.20] | 0.243 | 1.004 | 1346 | 2146 |
| Mean pressure experienced | 3.30* | 1.70 | [ 1.21, 9.23] | 0.990 | [0.84, 1.20] | 0.019 | 1.004 | 1610 | 2466 |
| Mean pressure utilized (partner’s view) | 1.23 | 0.70 | [ 0.41, 3.54] | 0.637 | [0.84, 1.20] | 0.226 | 1.003 | 1666 | 2712 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.78 | 0.19 | [0.46, 1.22] | NA | NA | NA | 1.001 | 1742 | 2619 |
| sd(Daily persuasion experienced) | 0.16 | 0.12 | [0.01, 0.44] | NA | NA | NA | 1.008 | 749 | 1628 |
| sd(Daily persuasion utilized (partner’s view)) | 0.18 | 0.13 | [0.01, 0.45] | NA | NA | NA | 1.002 | 1314 | 1689 |
| sd(Daily pressure experienced) | 0.47 | 0.22 | [0.07, 1.02] | NA | NA | NA | 1.007 | 1083 | 919 |
| sd(Daily pressure utilized (partner’s view)) | 0.31 | 0.31 | [0.02, 1.37] | NA | NA | NA | 1.002 | 1028 | 1907 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| disc | 1.00 | 0.00 | [1.00, 1.00] | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(reactance_ordinal),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
## = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.83),
## b_pressure_self_cb and b_persuasion_self_cb (r = 0.79),
## b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.84). This might
## lead to inappropriate results. See 'Details' in '?rope'.
conditional_spaghetti(
reactance_ordinal,
effects = c('persuasion_self_cw', 'pressure_self_cw')
, group_var = 'coupleID'
#, n_groups = 15
, plot_full_range = T
)\(persuasion_self_cw
<img
src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png"
width="2400" />\)pressure_self_cw
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}
formula <- bf(
is_reactance ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_nopushing_", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 31 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 756 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -361.7 15.6
## p_loo 63.7 4.9
## looic 723.3 31.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.5]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 725 95.9% 78
## (0.7, 1] (bad) 30 4.0% <NA>
## (1, Inf) (very bad) 1 0.1% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.002645503
summarize_brms(
is_reactance,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 0.32*** | 0.08 | [0.18, 0.53] | 1.000 | [0.83, 1.20] | 0.000 | 1.002 | 2308 | 2705 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.86 | 0.08 | [0.71, 1.01] | 0.964 | [0.83, 1.20] | 0.612 | 1.000 | 2887 | 2614 |
| Daily persuasion utilized (partner’s view) | 1.10 | 0.14 | [0.86, 1.45] | 0.777 | [0.83, 1.20] | 0.729 | 1.001 | 2245 | 2471 |
| Daily pressure experienced | 2.09* | 0.58 | [1.19, 4.20] | 0.991 | [0.83, 1.20] | 0.024 | 1.001 | 1589 | 1427 |
| Daily pressure utilized (partner’s view) | 1.34 | 0.53 | [0.57, 3.87] | 0.780 | [0.83, 1.20] | 0.276 | 1.001 | 1658 | 1539 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 1.86 | 0.67 | [0.89, 3.88] | 0.955 | [0.83, 1.20] | 0.092 | 1.001 | 4640 | 2592 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | 1.86 | 1.07 | [0.58, 5.93] | 0.859 | [0.83, 1.20] | 0.144 | 1.005 | 1143 | 2107 |
| Mean persuasion utilized (partner’s view) | 1.03 | 0.66 | [0.31, 3.66] | 0.520 | [0.83, 1.20] | 0.233 | 1.004 | 1265 | 2190 |
| Mean pressure experienced | 17.76*** | 16.40 | [3.37, 123.01] | 1.000 | [0.83, 1.20] | 0.000 | 1.002 | 1423 | 2478 |
| Mean pressure utilized (partner’s view) | 1.21 | 1.21 | [0.17, 7.77] | 0.574 | [0.83, 1.20] | 0.136 | 1.001 | 1712 | 2485 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 1.21 | 0.25 | [0.80, 1.78] | NA | NA | NA | 1.001 | 1396 | 1835 |
| sd(Daily persuasion experienced) | 0.19 | 0.13 | [0.01, 0.48] | NA | NA | NA | 1.004 | 958 | 1649 |
| sd(Daily persuasion utilized (partner’s view)) | 0.43 | 0.17 | [0.12, 0.83] | NA | NA | NA | 1.003 | 1279 | 1064 |
| sd(Daily pressure experienced) | 0.79 | 0.53 | [0.07, 2.12] | NA | NA | NA | 1.004 | 481 | 715 |
| sd(Daily pressure utilized (partner’s view)) | 0.74 | 0.64 | [0.04, 2.71] | NA | NA | NA | 1.002 | 1030 | 1680 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(is_reactance),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
process_model_component <- function(obj) {
# Convert to string, modify, and evaluate
name <- deparse(substitute(obj))
if (report_hurdle) name <- paste0(name, '_hu')
if (report_ordinal) name <- paste0(name, '_ordinal')
return(get(name, envir = parent.frame()))
}
all_models <- report_side_by_side(
pa_sub,
pa_obj_log,
mood_gauss,
reactance_ordinal,
is_reactance,
stats_to_report = c('CI', 'pd'),
model_rows_random = process_model_component(model_rows_random),
model_rows_fixed = process_model_component(model_rows_fixed),
model_rownames_random = process_model_component(model_rownames_random),
model_rownames_fixed = process_model_component(model_rownames_fixed)
) [1] “pa_sub”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.
[1] “pa_obj_log”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.
[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”
# pretty printing
summary_all_models <- all_models %>%
print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
add_header_above(
c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),
"Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5),
"Mood Gaussian" = (length(all_models) / 5),
"Reactance Ordinal" = (length(all_models) / 5),
"Reactance Dichotome" = (length(all_models) / 5))
)
export_xlsx(
summary_all_models,
rows_to_pack = process_model_component(rows_to_pack),
file.path("Output", "AllModelsNoPushing.xlsx"),
merge_option = 'header',
simplify_2nd_row = TRUE,
line_above_rows = c(1,2),
line_below_rows = c(-1)
)
print(summary_all_models)| exp(Est.) pa_sub | 95% CI pa_sub | pd pa_sub | exp(Est.) pa_obj_log | 95% CI pa_obj_log | pd pa_obj_log | Est. mood_gauss | 95% CI mood_gauss | pd mood_gauss | OR reactance_ordinal | 95% CI reactance_ordinal | pd reactance_ordinal | OR is_reactance | 95% CI is_reactance | pd is_reactance | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Intercept | 48.72*** | [42.92, 55.10] | 1.000 | 116.98*** | [105.53, 130.15] | 1.000 | 3.70*** | [ 3.49, 3.90] | 1.000 | NA | NA | NA | 0.32*** | [0.18, 0.53] | 1.000 |
| Hurdle Intercept | 0.82 | [ 0.59, 1.13] | 0.882 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Within-Person Effects | |||||||||||||||
| Daily persuasion experienced | 1.03 | [ 0.97, 1.09] | 0.853 | 1.04* | [ 1.01, 1.07] | 0.989 | 0.01 | [-0.03, 0.05] | 0.663 | 0.85* | [ 0.72, 0.99] | 0.979 | 0.86 | [0.71, 1.01] | 0.964 |
| Daily persuasion utilized (partner’s view) | 1.03 | [ 0.98, 1.08] | 0.877 | 1.02 | [ 0.99, 1.06] | 0.936 | 0.04 | [-0.01, 0.08] | 0.938 | 1.01 | [ 0.84, 1.20] | 0.548 | 1.10 | [0.86, 1.45] | 0.777 |
| Daily pressure experienced | 0.90* | [ 0.81, 0.99] | 0.983 | 0.95 | [ 0.89, 1.02] | 0.935 | -0.03 | [-0.14, 0.07] | 0.743 | 1.98** | [ 1.31, 2.76] | 0.999 | 2.09* | [1.19, 4.20] | 0.991 |
| Daily pressure utilized (partner’s view) | 0.94 | [ 0.86, 1.02] | 0.928 | 0.99 | [ 0.93, 1.06] | 0.627 | 0.01 | [-0.10, 0.11] | 0.559 | 1.15 | [ 0.72, 1.85] | 0.749 | 1.34 | [0.57, 3.87] | 0.780 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.98 | [ 0.87, 1.11] | 0.623 | 0.97 | [ 0.91, 1.04] | 0.774 | 0.25*** | [ 0.15, 0.36] | 1.000 | 1.63 | [ 0.87, 3.10] | 0.929 | 1.86 | [0.89, 3.88] | 0.955 |
| Daily weartime | NA | NA | NA | 1.00*** | [ 1.00, 1.00] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||||||||
| Mean persuasion experienced | 1.11 | [ 0.85, 1.47] | 0.775 | 1.11 | [ 0.85, 1.42] | 0.788 | 0.42 | [-0.10, 0.91] | 0.942 | 1.08 | [ 0.42, 2.76] | 0.567 | 1.86 | [0.58, 5.93] | 0.859 |
| Mean persuasion utilized (partner’s view) | 1.09 | [ 0.83, 1.45] | 0.737 | 1.04 | [ 0.79, 1.34] | 0.611 | 0.33 | [-0.18, 0.82] | 0.902 | 0.79 | [ 0.29, 2.26] | 0.665 | 1.03 | [0.31, 3.66] | 0.520 |
| Mean pressure experienced | 1.12 | [ 0.81, 1.55] | 0.752 | 0.89 | [ 0.68, 1.20] | 0.782 | -0.36 | [-0.91, 0.19] | 0.903 | 3.30* | [ 1.21, 9.23] | 0.990 | 17.76*** | [3.37, 123.01] | 1.000 |
| Mean pressure utilized (partner’s view) | 0.91 | [ 0.64, 1.28] | 0.708 | 1.02 | [ 0.78, 1.37] | 0.578 | -0.29 | [-0.82, 0.27] | 0.851 | 1.23 | [ 0.41, 3.54] | 0.637 | 1.21 | [0.17, 7.77] | 0.574 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | 1.00 | [ 1.00, 1.00] | 0.912 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||||||||
| Hu Daily persuasion experienced | 1.61*** | [ 1.43, 1.84] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | 1.43*** | [ 1.28, 1.61] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | 1.01 | [ 0.76, 1.39] | 0.525 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | 1.74*** | [ 1.26, 2.83] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.92 | [ 0.71, 1.20] | 0.730 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||||||||
| Hu Mean persuasion experienced | 1.56 | [ 0.76, 3.18] | 0.889 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | 1.55 | [ 0.76, 3.19] | 0.881 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | 0.33** | [ 0.15, 0.77] | 0.996 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | 0.58 | [ 0.25, 1.33] | 0.911 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||||||||
| sd(Intercept) | 0.32 | [0.24, 0.42] | NA | 0.30 | [0.23, 0.39] | NA | 0.61 | [0.48, 0.79] | NA | 0.78 | [0.46, 1.22] | NA | 1.21 | [0.80, 1.78] | NA |
| sd(Hurdle Intercept) | 0.89 | [0.70, 1.18] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion experienced) | 0.12 | [0.08, 0.17] | NA | 0.05 | [0.03, 0.08] | NA | 0.05 | [0.00, 0.11] | NA | 0.16 | [0.01, 0.44] | NA | 0.19 | [0.01, 0.48] | NA |
| sd(Daily persuasion utilized (partner’s view)) | 0.09 | [0.05, 0.14] | NA | 0.06 | [0.03, 0.09] | NA | 0.08 | [0.03, 0.14] | NA | 0.18 | [0.01, 0.45] | NA | 0.43 | [0.12, 0.83] | NA |
| sd(Daily pressure experienced) | 0.08 | [0.00, 0.25] | NA | 0.05 | [0.00, 0.14] | NA | 0.06 | [0.00, 0.23] | NA | 0.47 | [0.07, 1.02] | NA | 0.79 | [0.07, 2.12] | NA |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | [0.00, 0.19] | NA | 0.04 | [0.00, 0.12] | NA | 0.07 | [0.00, 0.25] | NA | 0.31 | [0.02, 1.37] | NA | 0.74 | [0.04, 2.71] | NA |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion experienced) | 0.21 | [0.08, 0.37] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | 0.19 | [0.05, 0.34] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | 0.23 | [0.01, 0.80] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | 0.31 | [0.02, 1.01] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||||||||
| sigma | 0.69 | [0.67, 0.71] | NA | 0.58 | [0.56, 0.59] | NA | 0.96 | [0.94, 0.98] | NA | NA | NA | NA | NA | NA | NA |
formula <- bf(
pa_sub ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(1 + pushing_self_cw + pushing_partner_cw | coupleID),
hu = ~ pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(1 + pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
, brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_onlypushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
## Start sampling
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 1342 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -5692.9 78.2
## p_loo 97.8 4.1
## looic 11385.8 156.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 1341 99.9% 633
## (0.7, 1] (bad) 1 0.1% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 4, observations = 1342, p-value = 0.26
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.003371833
## sample estimates:
## outlier frequency (expected: 0.00149776453055142 )
## 0.002980626
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
pa_sub,
stats_to_report = c('pd','ROPE'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub_continuous <- summarize_brms(
pa_sub,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_continuous,
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous
columns_to_replace <- c("ROPE", "inside ROPE")
summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <-
summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]
# Print the updated dataframe
summary_pa_sub %>%
print_df(rows_to_pack = rows_to_pack_hu)| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 57.28*** | 4.53 | [49.10, 66.97] | 1.000 | [0.93, 1.08] | 0.000 | 1.002 | 1327 | 2040 |
| Hurdle Intercept | 3.79*** | 0.98 | [ 2.27, 6.26] | 1.000 | [0.84, 1.20] | 0.000 | 1.005 | 1088 | 2463 |
| Conditional Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.97 | 0.03 | [ 0.91, 1.03] | 0.856 | [0.93, 1.08] | 0.931 | 1.002 | 4427 | 2662 |
| Daily pushing utilized (partner’s view) | 0.97 | 0.03 | [ 0.90, 1.03] | 0.847 | [0.93, 1.08] | 0.907 | 1.001 | 3679 | 2873 |
| Day | 0.91 | 0.07 | [ 0.79, 1.06] | 0.894 | [0.93, 1.08] | 0.398 | 1.000 | 8744 | 2600 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 1.32* | 0.18 | [ 1.00, 1.71] | 0.977 | [0.93, 1.08] | 0.060 | 1.005 | 1217 | 2176 |
| Mean pushing utilized (partner’s view) | 1.18 | 0.15 | [ 0.90, 1.53] | 0.883 | [0.93, 1.08] | 0.222 | 1.005 | 1265 | 2076 |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||
| Hu Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | 1.55** | 0.23 | [ 1.17, 2.24] | 0.998 | [0.84, 1.20] | 0.035 | 1.002 | 3127 | 2531 |
| Hu Daily pushing utilized (partner’s view) | 1.84*** | 0.30 | [ 1.38, 2.68] | 1.000 | [0.84, 1.20] | 0.002 | 1.001 | 2766 | 2989 |
| Hu Day | 0.61 | 0.15 | [ 0.37, 1.00] | 0.974 | [0.84, 1.20] | 0.100 | 1.003 | 7632 | 2768 |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||
| Hu Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | 0.46** | 0.14 | [ 0.26, 0.81] | 0.997 | [0.84, 1.20] | 0.018 | 1.001 | 2675 | 2824 |
| Hu Mean pushing utilized (partner’s view) | 0.62 | 0.18 | [ 0.33, 1.08] | 0.952 | [0.84, 1.20] | 0.136 | 1.001 | 2568 | 2699 |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.37 | 0.06 | [0.29, 0.50] | NA | NA | NA | 1.002 | 1283 | 2370 |
| sd(Hurdle Intercept) | 1.25 | 0.19 | [0.93, 1.74] | NA | NA | NA | 1.002 | 1448 | 2426 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.06 | 0.04 | [0.00, 0.16] | NA | NA | NA | 1.001 | 870 | 1486 |
| sd(Daily pushing utilized (partner’s view)) | 0.09 | 0.04 | [0.01, 0.18] | NA | NA | NA | 1.004 | 977 | 1142 |
| sd(Hu Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | 0.46 | 0.24 | [0.07, 1.02] | NA | NA | NA | 1.000 | 947 | 1689 |
| sd(Hu Daily pushing utilized (partner’s view)) | 0.51 | 0.24 | [0.10, 1.13] | NA | NA | NA | 1.001 | 1088 | 1582 |
| Additional Parameters | |||||||||
| sigma | 0.67 | 0.02 | [0.64, 0.71] | NA | NA | NA | 1.002 | 6374 | 2773 |
# Plot continuous part of model
variable <- c(
'(Intercept)',
'b_pushing_self_cw',
'b_pushing_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()plot(
bayestestR::rope(
pa_sub,
parameter = variable,
range = rope_range_continuous,
verbose = F,
ci = 1
)
) + theme_bw()# Hurdle part of the model
variable <- c(
'b_hu_pushing_self_cw',
'b_hu_pushing_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
verbose = FALSE
) + theme_bw()## Possible multicollinearity between b_pushing_partner_cb and
## b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
## See 'Details' in '?rope'.
conds_eff <- conditional_spaghetti(
pa_sub,
effects = c(
'pushing_self_cw',
'pushing_partner_cw'
),
x_label = c(
'Received Plan-Related Pushing',
'Exerted Plan-Related Pushing'
),
group_var = 'coupleID',
plot_full_range = TRUE,
y_limits = c(0, 100),
y_label = "Same-Day MVPA",
y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
, filter_quantiles = .9995
, font_family = 'Candara'
)## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
$pushing_self_cw
## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0172
$pushing_partner_cw
## Warning: Removed 130 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 180 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.022
Note. This graphic illustrates the relationship between
social control and moderate to vigorous physical activity (MVPA) using a
Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered
within individuals to examine how deviations from their average social
control relate to same-day MVPA. Shaded areas indicate credible
intervals, thick lines show fixed effects, and thin lines represent
random effects, highlighting variability across couples. The plots
display the probability of being active, expected minutes of MVPA when
active, and combined predicted MVPA. The bottom density plot visualizes
the posterior distributions of slope estimates, transformed to represent
multiplicative changes in odds ratios (hurdle component) or expected
values. Medians and 95% credible intervals (2.5th and 97.5th
percentiles) are shown. Effects are significant, when the 95% credible
interval does not overlap 1.
## [1] 5.75 971.25
formula <- bf(
pa_obj ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day + weartime_self_cw + weartime_self_cb +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "Intercept")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_onlypushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## Computed from 4000 by 1214 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -6914.2 39.4
## p_loo 58.6 4.3
## looic 13828.4 78.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 8, observations = 1214, p-value = 0.16
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.0008237232 0.0074135091
## sample estimates:
## outlier frequency (expected: 0.00339373970345964 )
## 0.006589786
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summarize_brms(
pa_obj_log,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_log,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.
| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 124.24*** | 7.52 | [109.69, 139.58] | 1.000 | [0.94, 1.06] | 0.000 | 1.002 | 967 | 1634 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.03 | 0.03 | [ 0.97, 1.08] | 0.807 | [0.94, 1.06] | 0.914 | 1.001 | 3221 | 2733 |
| Daily pushing utilized (partner’s view) | 1.02 | 0.02 | [ 0.97, 1.06] | 0.797 | [0.94, 1.06] | 0.979 | 1.000 | 5065 | 2748 |
| Day | 0.97 | 0.05 | [ 0.87, 1.08] | 0.718 | [0.94, 1.06] | 0.654 | 1.000 | 7350 | 3336 |
| Daily weartime | 1.00** | 0.00 | [ 1.00, 1.00] | 0.998 | [0.94, 1.06] | 1.000 | 1.001 | 4680 | 2898 |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 0.97 | 0.06 | [ 0.86, 1.10] | 0.680 | [0.94, 1.06] | 0.635 | 1.001 | 2223 | 2645 |
| Mean pushing utilized (partner’s view) | 1.04 | 0.07 | [ 0.92, 1.18] | 0.739 | [0.94, 1.06] | 0.578 | 1.002 | 2058 | 2692 |
| Mean weartime | 1.00 | 0.00 | [ 1.00, 1.00] | 0.812 | [0.94, 1.06] | 1.000 | 1.000 | 3879 | 3246 |
| Random Effects | |||||||||
| sd(Intercept) | 0.31 | 0.04 | [0.24, 0.42] | NA | NA | NA | 1.005 | 890 | 1790 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.09 | 0.04 | [0.02, 0.17] | NA | NA | NA | 1.004 | 887 | 821 |
| sd(Daily pushing utilized (partner’s view)) | 0.03 | 0.03 | [0.00, 0.10] | NA | NA | NA | 1.001 | 1398 | 1731 |
| Additional Parameters | |||||||||
| sigma | 0.54 | 0.01 | [0.52, 0.57] | NA | NA | NA | 1.001 | 6040 | 3037 |
plot(
bayestestR::p_direction(pa_obj_log),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## [1] 0 5
formula <- bf(
aff ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss_onlypushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## Computed from 4000 by 1342 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -1798.8 39.9
## p_loo 50.3 3.0
## looic 3597.5 79.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.8]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 13, observations = 1342, p-value =
## 4.837e-06
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.005167734 0.016508158
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.009687034
summarize_brms(
mood_gauss,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)| Est. | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 3.78*** | 0.12 | [ 3.56, 4.00] | 1.000 | [-0.11, 0.11] | 0.000 | 1.013 | 560 | 1111 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.02 | 0.04 | [-0.05, 0.08] | 0.670 | [-0.11, 0.11] | 0.995 | 1.002 | 2823 | 2591 |
| Daily pushing utilized (partner’s view) | 0.09* | 0.04 | [ 0.01, 0.17] | 0.984 | [-0.11, 0.11] | 0.715 | 1.002 | 1597 | 2207 |
| Day | 0.28*** | 0.09 | [ 0.11, 0.46] | 1.000 | [-0.11, 0.11] | 0.023 | 1.001 | 3871 | 3088 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 0.12 | 0.11 | [-0.08, 0.32] | 0.880 | [-0.11, 0.11] | 0.433 | 1.003 | 1269 | 1954 |
| Mean pushing utilized (partner’s view) | 0.08 | 0.11 | [-0.12, 0.30] | 0.786 | [-0.11, 0.11] | 0.554 | 1.002 | 1216 | 1958 |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.61 | 0.08 | [0.48, 0.79] | NA | NA | NA | 1.011 | 648 | 1448 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.07 | 0.04 | [0.00, 0.16] | NA | NA | NA | 1.000 | 1409 | 1126 |
| sd(Daily pushing utilized (partner’s view)) | 0.14 | 0.04 | [0.06, 0.24] | NA | NA | NA | 1.002 | 1329 | 828 |
| Additional Parameters | |||||||||
| sigma | 0.91 | 0.02 | [0.87, 0.94] | NA | NA | NA | 1.001 | 3887 | 2536 |
plot(
bayestestR::p_direction(mood_gauss),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
conditional_spaghetti(
mood_gauss,
effects = c('pushing_partner_cw'),
group_var = 'coupleID',
plot_full_range = TRUE
)$pushing_partner_cw
## [1] 0 5
df_double$reactance_ordinal <- factor(df_double$reactance,
levels = 0:5,
ordered = TRUE)
formula <- bf(
reactance_ordinal ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = cumulative() # HURDLE_CUMULATIVE
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance_ordinal <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::cumulative(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777
, file = file.path("models_cache_brms", paste0("reactance_ordinal_onlypushing", suffix))
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
##
## Computed from 4000 by 403 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -363.6 23.5
## p_loo 33.6 2.9
## looic 727.1 47.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.6]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## DHARMa bootstrapped outlier test
##
## data: model.check
## outliers at both margin(s) = 0, observations = 403, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.00248139
## sample estimates:
## outlier frequency (expected: 0.000297766749379653 )
## 0
summarize_brms(
reactance_ordinal,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_ordinal,
model_rows_random = model_rows_random_ordinal,
model_rownames_fixed = model_rownames_fixed_ordinal,
model_rownames_random = model_rownames_random_ordinal,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack_ordinal)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercepts | |||||||||
| Intercept | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Intercept[1] | 5.44*** | 1.57 | [ 3.06, 10.06] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3521 | 2882 |
| Intercept[2] | 9.89*** | 3.14 | [ 5.45, 18.74] | 1.000 | [0.84, 1.20] | 0.000 | 1.002 | 3598 | 2951 |
| Intercept[3] | 23.10*** | 7.89 | [ 12.38, 46.79] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3933 | 2899 |
| Intercept[4] | 81.88*** | 35.09 | [ 38.53, 199.63] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 4569 | 3054 |
| Intercept[5] | 771.25*** | 625.85 | [202.52, 4855.18] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 5246 | 2545 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.31* | 0.14 | [ 1.04, 1.62] | 0.987 | [0.84, 1.20] | 0.214 | 1.001 | 3641 | 2884 |
| Daily pushing utilized (partner’s view) | 0.96 | 0.13 | [ 0.75, 1.25] | 0.615 | [0.84, 1.20] | 0.810 | 1.002 | 4044 | 2794 |
| Day | 1.60 | 0.73 | [ 0.68, 3.95] | 0.854 | [0.84, 1.20] | 0.182 | 1.002 | 4822 | 3245 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 2.32** | 0.65 | [ 1.35, 3.99] | 1.000 | [0.84, 1.20] | 0.005 | 1.001 | 3394 | 2987 |
| Mean pushing utilized (partner’s view) | 1.25 | 0.35 | [ 0.70, 2.16] | 0.783 | [0.84, 1.20] | 0.365 | 1.001 | 3738 | 2562 |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.85 | 0.21 | [0.49, 1.34] | NA | NA | NA | 1.001 | 1569 | 2204 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.17 | 0.14 | [0.01, 0.50] | NA | NA | NA | 1.004 | 1241 | 1651 |
| sd(Daily pushing utilized (partner’s view)) | 0.16 | 0.14 | [0.01, 0.57] | NA | NA | NA | 1.003 | 1626 | 1714 |
| Additional Parameters | |||||||||
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| disc | 1.00 | 0.00 | [1.00, 1.00] | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(reactance_ordinal),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
## = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
## to inappropriate results. See 'Details' in '?rope'.
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}
formula <- bf(
is_reactance ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance_onlypushing", suffix))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.
##
## Divergences:
## 0 of 4000 iterations ended with a divergence.
##
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
##
## Energy:
## E-BFMI indicated no pathological behavior.
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
## Using 10 posterior draws for ppc type 'dens_overlay' by default.
## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
##
## Computed from 4000 by 403 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -216.1 11.0
## p_loo 39.3 3.0
## looic 432.3 22.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 2.0]).
##
## Pareto k diagnostic values:
## Count Pct. Min. ESS
## (-Inf, 0.7] (good) 400 99.3% 458
## (0.7, 1] (bad) 3 0.7% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: model.check
## outliers at both margin(s) = 1, observations = 403, p-value = 0.5534
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 6.282137e-05 1.374725e-02
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.00248139
summarize_brms(
is_reactance,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 0.40** | 0.11 | [0.23, 0.68] | 1.000 | [0.83, 1.20] | 0.002 | 1.000 | 3095 | 2830 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.37* | 0.18 | [1.06, 1.87] | 0.993 | [0.83, 1.20] | 0.154 | 1.000 | 3238 | 2410 |
| Daily pushing utilized (partner’s view) | 1.16 | 0.26 | [0.80, 2.12] | 0.769 | [0.83, 1.20] | 0.509 | 1.000 | 2310 | 2170 |
| Day | 1.78 | 0.89 | [0.69, 4.81] | 0.874 | [0.83, 1.20] | 0.160 | 1.000 | 6984 | 2927 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 5.39*** | 2.58 | [2.38, 14.94] | 1.000 | [0.83, 1.20] | 0.000 | 1.001 | 2293 | 2628 |
| Mean pushing utilized (partner’s view) | 2.37 | 1.15 | [0.97, 6.82] | 0.972 | [0.83, 1.20] | 0.057 | 1.001 | 2556 | 2640 |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 1.52 | 0.31 | [1.03, 2.24] | NA | NA | NA | 1.001 | 1431 | 2494 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.18 | 0.15 | [0.01, 0.56] | NA | NA | NA | 1.003 | 1256 | 1876 |
| sd(Daily pushing utilized (partner’s view)) | 0.42 | 0.27 | [0.03, 1.15] | NA | NA | NA | 1.002 | 1292 | 1879 |
| Additional Parameters | |||||||||
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(is_reactance),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
conditional_spaghetti(
is_reactance,
effects = c('pushing_self_cw'),
group_var = 'coupleID',
plot_full_range = TRUE
)$pushing_self_cw
process_model_component <- function(obj) {
# Convert to string, modify, and evaluate
name <- deparse(substitute(obj))
if (report_hurdle) name <- paste0(name, '_hu')
if (report_ordinal) name <- paste0(name, '_ordinal')
return(get(name, envir = parent.frame()))
}
all_models <- report_side_by_side(
pa_sub,
pa_obj_log,
mood_gauss,
reactance_ordinal,
is_reactance,
stats_to_report = c('CI', 'pd'),
model_rows_random = process_model_component(model_rows_random),
model_rows_fixed = process_model_component(model_rows_fixed),
model_rownames_random = process_model_component(model_rownames_random),
model_rownames_fixed = process_model_component(model_rownames_fixed)
) [1] “pa_sub”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.
[1] “pa_obj_log”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.
[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”
# pretty printing
summary_all_models <- all_models %>%
print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
add_header_above(
c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),
"Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5),
"Mood Gaussian" = (length(all_models) / 5),
"Reactance Ordinal" = (length(all_models) / 5),
"Reactance Dichotome" = (length(all_models) / 5))
)
export_xlsx(
summary_all_models,
rows_to_pack = process_model_component(rows_to_pack),
file.path("Output", "AllModelsOnlyPushing.xlsx"),
merge_option = 'header',
simplify_2nd_row = TRUE,
line_above_rows = c(1,2),
line_below_rows = c(-1)
)
print(summary_all_models)| exp(Est.) pa_sub | 95% CI pa_sub | pd pa_sub | exp(Est.) pa_obj_log | 95% CI pa_obj_log | pd pa_obj_log | Est. mood_gauss | 95% CI mood_gauss | pd mood_gauss | OR reactance_ordinal | 95% CI reactance_ordinal | pd reactance_ordinal | OR is_reactance | 95% CI is_reactance | pd is_reactance | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Intercept | 57.28*** | [49.10, 66.97] | 1.000 | 124.24*** | [109.69, 139.58] | 1.000 | 3.78*** | [ 3.56, 4.00] | 1.000 | NA | NA | NA | 0.40** | [0.23, 0.68] | 1.000 |
| Hurdle Intercept | 3.79*** | [ 2.27, 6.26] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Within-Person Effects | |||||||||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.97 | [ 0.91, 1.03] | 0.856 | 1.03 | [ 0.97, 1.08] | 0.807 | 0.02 | [-0.05, 0.08] | 0.670 | 1.31* | [ 1.04, 1.62] | 0.987 | 1.37* | [1.06, 1.87] | 0.993 |
| Daily pushing utilized (partner’s view) | 0.97 | [ 0.90, 1.03] | 0.847 | 1.02 | [ 0.97, 1.06] | 0.797 | 0.09* | [ 0.01, 0.17] | 0.984 | 0.96 | [ 0.75, 1.25] | 0.615 | 1.16 | [0.80, 2.12] | 0.769 |
| Day | 0.91 | [ 0.79, 1.06] | 0.894 | 0.97 | [ 0.87, 1.08] | 0.718 | 0.28*** | [ 0.11, 0.46] | 1.000 | 1.60 | [ 0.68, 3.95] | 0.854 | 1.78 | [0.69, 4.81] | 0.874 |
| Daily weartime | NA | NA | NA | 1.00** | [ 1.00, 1.00] | 0.998 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 1.32* | [ 1.00, 1.71] | 0.977 | 0.97 | [ 0.86, 1.10] | 0.680 | 0.12 | [-0.08, 0.32] | 0.880 | 2.32** | [ 1.35, 3.99] | 1.000 | 5.39*** | [2.38, 14.94] | 1.000 |
| Mean pushing utilized (partner’s view) | 1.18 | [ 0.90, 1.53] | 0.883 | 1.04 | [ 0.92, 1.18] | 0.739 | 0.08 | [-0.12, 0.30] | 0.786 | 1.25 | [ 0.70, 2.16] | 0.783 | 2.37 | [0.97, 6.82] | 0.972 |
| Mean weartime | NA | NA | NA | 1.00 | [ 1.00, 1.00] | 0.812 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||||||||
| Hu Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | 1.55** | [ 1.17, 2.24] | 0.998 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | 1.84*** | [ 1.38, 2.68] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.61 | [ 0.37, 1.00] | 0.974 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||||||||
| Hu Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | 0.46** | [ 0.26, 0.81] | 0.997 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | 0.62 | [ 0.33, 1.08] | 0.952 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||||||||
| sd(Intercept) | 0.37 | [0.29, 0.50] | NA | 0.31 | [0.24, 0.42] | NA | 0.61 | [0.48, 0.79] | NA | 0.85 | [0.49, 1.34] | NA | 1.52 | [1.03, 2.24] | NA |
| sd(Hurdle Intercept) | 1.25 | [0.93, 1.74] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.06 | [0.00, 0.16] | NA | 0.09 | [0.02, 0.17] | NA | 0.07 | [0.00, 0.16] | NA | 0.17 | [0.01, 0.50] | NA | 0.18 | [0.01, 0.56] | NA |
| sd(Daily pushing utilized (partner’s view)) | 0.09 | [0.01, 0.18] | NA | 0.03 | [0.00, 0.10] | NA | 0.14 | [0.06, 0.24] | NA | 0.16 | [0.01, 0.57] | NA | 0.42 | [0.03, 1.15] | NA |
| sd(Hu Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | 0.46 | [0.07, 1.02] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | 0.51 | [0.10, 1.13] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||||||||
| sigma | 0.67 | [0.64, 0.71] | NA | 0.54 | [0.52, 0.57] | NA | 0.91 | [0.87, 0.94] | NA | NA | NA | NA | NA | NA | NA |
Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)